Critical Boolean networks with scale-free in-degree distribution 
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We investigate analytically and numerically the dynamical properties of critical Boolean networks 
with power-law in-degree distributions. When the exponent of the in-degree distribution is larger 
than 3, we obtain results equivalent to those obtained for networks with fixed in-degree, e.g., the 
number of the non-frozen nodes scales as N 2 ^ 3 with the system size N. When the exponent of the 
distribution is between 2 and 3, the number of the non-frozen nodes increases as N x , with x being 
between and 2/3 and depending on the exponent and on the cutoff of the in-degree distribution. 
These and ensuing results explain various findings obtained earlier by computer simulations. 
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Complex dynamical systems, where a large number of 
units interact in a non-trivial way, are often modelled as 
networks. The units from which these networks are built 
can show various types of intrinsic dynamics, including 
oscillations. Whenever the dynamics can be reduced to 
only two possible states per node, a Boolean network is 
obtained. Stuart Kauffman was the first to use Ran- 
dom Boolean networks (RBNs) to model the dynamics 
of genetic and protein networks |l|, l2| . Although Boolean 
models represent a strong simplification of the far more 
complex reality, there exist several examples where the 
modelling of a cellular network by Boolean variables cap- 
tures correctly the essential dynamics of the system 0, [j] . 

RBNs are directed graphs where each node i has a 
Boolean value <7i £ {0,1} and an update function 
which determines the new value in the next time step 
as function of the state of those nodes that have a link 
to node i. Links and functions are assigned at random, 
given certain constraints concerning the number of in- 
puts per node or the set of functions. The update can be 
performed in different ways, we consider here the usual 
case of synchronous update. After some time, the dy- 
namics reaches an attractor, i.e. a periodic sequence of 
states. Depending on the parameters of the network, the 
dynamics is either in the frozen phase or in the chaotic 
phase or at the critical point between the two. In the 
frozen phase, all apart from a small number of nodes as- 
sume a constant value on the attractors, i.e., they are 
frozen. When the state of a node is changed, on average 
less than one node will be changed in the next time step, 
and the size of a perturbation decreases with time. In the 
chaotic phase, a nonvanishing proportion of nodes keep 
changing their state even after a long time. The size of 
a perturbation increases with time, since a change in the 
state of one node will lead on an average to a change of 
the state of more than one nodes in the next time step. 
Most studies of RBNs have focused on the critical point, 
which is at the boundary between these two phases, and 
where a perturbation of one node propagates on an av- 
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erage to one other node. These studies deal mainly with 
the (mean) number and size of attractors, motivated by 
Kauffman's original claim that biological networks are 
poised at the critical point, and that attractors can be 
equated with cell types. Despite of the long time since the 
introduction of the model, a full analytical understanding 
of critical RBNs was obtained only recently [1, H, 0] • 

A key concept at understanding the dynamics of crit- 
ical RBNs is the classification of the nodes according to 
their dynamical behavior on attractors into frozen, non- 
frozen and relevant nodes 0]. Relevant nodes are those 
nodes that determine the attractors, while the other non- 
frozen nodes are slaved to the dynamics of the relevant 
nodes; changing their state does not change the attrac- 
tor. A stochastic process that gradually determines the 
frozen core starting from the nodes that have a constant 
function, was used in (i. ;9j to prove that the number 
of nonfrozen nodes in critical RBNs scales as N 2 ^ 3 , and 
the number of relevant nodes as TV 1 / 3 , with TV being the 
number of nodes in the network. In the limit of large net- 
work size, scaling functions for the number of non-frozen 
and relevant nodes were calculated analytically. 

All studies mentioned so far assign k inputs to each 
node, while the number of outputs is Poisson distributed, 
since incoming links are connected at random to a node 
where they originate. However, biological networks are 
known to have a broad degree distribution, which is of- 
ten well described by a power law (see [Io| and references 
therein). For this reason, several recent studies were de- 
voted to Boolean dynamics on scale-free networks. The 
majority of these studies uses a scale-free in-degree dis- 
tribution and a Poissonian out-degree distribution, but 
other implementations can also be found. Observations 
made in computer simulations are that attractors are 
shorter and frozen nodes are more numerous in critical 
scale-free networks compared to RBNs with a fixed num- 
ber of in put s, given the same total number of links and 
of nodes [ll|, UM > an( i that attractors are sensitive to per- 
turbations of highly connected nodes, but not of sparsely 
connected nodes [l2|, [l3[ . These and other [3, EH sim- 
ulation results are merely stated and are not embedded 
into an analytical framework. Analytical results obtained 
so far are limited to calculating the phase diagram using 
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the annealed approximation [TH, [l6|, [TtJ ; only the work 
by Lee and Rieger [l8[ goes further by calculating the 
asymptotic Hamming distance in the chaotic phase and 
extrapolating the results to the critical point by using a 
finite-size scaling ansatz in combination with the calcu- 
lation of the size distribution of perturbed clusters. 

In this paper, we will present an analytical calcula- 
tion for RBNs with scale-free input distributions at the 
critical point, obtaining scaling laws for the number of 
nonfrozen and relevant nodes. Our results, which are 
confirmed by a numerical evaluation, explain the above- 
mentioned findings of computer simulations, and convey 
a clear understanding of the properties of attractors in 
these systems. 

We consider critical networks that have an in-degree 
distribution P(k) that follows a power law, P(k) = A-k -1 
for k > 1. The normalization constant A depends on 
the minimum and the maximum in-degree. We fixed 
the minimum in-degree to 2; the maximum in-degree de- 
pends on the network size and the chosen implementa- 
tion of the model (see below) . We consider only the case 
7 > 2, where such a normalization is possible. In the 
case 2 < 7 < 3, the second moment of the degree distri- 
bution diverges, and it has been argued in [1 81 ] that this 
should change the dynamical properties. The out-degree 
distribution is Poissonian with a mean (k) — J^k kP(k), 
since the input connections are chosen at random from 
all nodes, just as for RBNs with fixed k. 

We investigated two ways of creating the input dis- 
tribution. First, we assigned to each node i a number 
ki of inputs that was drawn from the distribution P(k), 
not allowing values fc, larger than N or smaller than 2. 
The total number of links and the largest value of ki 
differ in this case between different networks. Second, 
we fixed the number of nodes with k inputs exactly at 
the value NP(k) (rounded to the nearest integer), which 
gives a distribution P(fc) that has a cutoff at k ~ N 1 ^ 1 . 
In part of the above-mentioned studies, networks with 
scale-free in-degree distributions were generated using a 
constraint that does not allow multiple connections be- 
tween the same nodes, or using a preferential-attachment 
algorithm, however, all these are known to create correla- 
tions between the degree of neighboring nodes , which 
in turn can affect the dynamics on these networks [20| . 
In order to avoid such complications, we connect the in- 
coming links at random to any node, without imposing 
any constraints. 

We also investigated several ways of assigning the 
Boolean functions to the nodes. First, we chose biased 
functions with a parameter p, assigning to each of the 
2 ki input configurations the output 1 with a probability 
p and the output with a probability 1 — p. The value 
of p was chosen such that the network is critical, i.e., 
that p = l/(k) [13]. The main results did not depend on 
whether we chose the exact mean (which can be differ- 
ent for each network), or the theoretical mean J^k kP(k). 
The second way of assigning the Boolean functions is to 
take only constant and reversible functions. There are 



two constant functions, which fix the value of a node 
to either or 1, irrespective of its input values. For 
each value of fc, there are 2 reversible functions, which 
are defined by the condition that changing the value of 
one input always changes the output. A node with a re- 
versible function becomes frozen only if all of its inputs 
are frozen. Such a network is critical if the total number 
of nodes equals the total number of inputs to nodes with 
reversible functions. Links to nodes with constant func- 
tions have no effect and can be omitted, so that the total 
number of links becomes identical to the total number of 
nodes. 
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FIG. 1: A sketch of one step of the stochastic process to 
determine the frozen core, (a) Nodes are placed in containers 
according to the number of inputs of which we do not yet 
already know for sure that they are frozen. We choose a 
node o from the container Co- (b) This node becomes an 
input to a node in container k with probability k/N . In this 
example, it becomes the input to 2 nodes, (c) The frozen links 
are removed, the two nodes are moved to container Ck-i, and 
node o is removed from the system. 

We adjusted the method proposed in [|| in order to 
determine the size of the frozen core of critical networks 
with scale- free input distributions. In the following, we 
describe this method for the case of only constant and 
reversible functions. The generalization to other cases is 
straightforward. The frozen core is determined starting 
from the nodes with constant functions and determining 
stepwise all those nodes that become frozen because all 
their inputs are frozen. The main idea of our method 
is to not specify the network in advance, but to choose 
the connections within the network while determining the 
frozen core. To this purpose, we place all N nodes of the 
network into "containers" Ck according to the number k 
of inputs. As mentioned above, inputs to nodes with con- 
stant functions are omitted, and these nodes are there- 
fore put into container Cq. The largest container index 
is fcmax = N or fc max ~ N 1 / 1 ', depending on the method 
chosen for creating the input distribution. The contents 
of the containers change with time, since we remove step- 
wise all those inputs of which we know that they come 
from a frozen node. The "time" we are defining here is 
not the real time for the dynamics of the system, but 
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it counts the steps of the stochastic process that we use 
to determine the frozen core. During one time step, we 
choose one node from the container Co and determine to 
which nodes this node is an input. Since the inputs are 
picked at random, the chosen node is connected to each 
input with probability 1/N. These inputs are removed, 
and the corresponding nodes moved from container Ck to 
container C/-_i (or to a lower container, when more than 
one connection is made to the same node) . At the end of 
the time step, we remove the chosen node from the sys- 
tem, and the number A of nodes in the system is reduced 
by 1. Thus, at each time t the number \Ck\ of nodes in 
container Ck is the number of nodes that have k inputs 
that have not yet become frozen during the process. In 
container Co are those nodes of which we know already 
that they are frozen, but we have not yet determined to 
which other nodes they are an input. We denote from 
now on the total number of nodes in the system by Ai„i, 
which is identical to N(t = 0). The number A; n ; — N(t) 
of nodes have been removed from the system. They are 
those nodes for which we have already determined that 
they are frozen and to which other nodes they are an 
input. 

The process ends when there are no nodes left in con- 
tainer Co, or when all nodes are in container Co- In the 
latter case, the entire network freezes, and the dynamics 
of the system runs to the same fixed point for all initial 
conditions. In the first case, there is a set of nonfrozen 
nodes. In order to determine the topology of the non- 
frozen part of the network, one can then fix the connec- 
tions that have not yet been determined by connecting 
the remaining inputs at random to the remaining nodes. 

Before showing the results of our computer simulations 
of this process obtained for an ensemble of many net- 
works, let us first perform an analytical calculation in 
order to predict the mean number of nodes remaining in 
the different containers at the end. We begin by evalu- 
ating the mean number of nodes in container Ck at the 
moment where only the fraction e = N/N- ln i nodes are 
left in the system. At this moment, container Ck con- 
tains all nodes that had initially I > k inputs, and where 
l — k inputs have already become frozen. The probability 
that an input has not yet become frozen is identical to 
e, since only the proportion e of nodes have not yet been 
removed, and since an input is connected to every node 
with the same probability. Since container C; contained 
initially oc Ni n il~ 7 nodes, we have 

k 

\C k \ oc N ini |] r\\l - e)^ fc Q • (1) 

For small e, nodes in container C& originated in contain- 
ers C with I ^> k, and we can therefore set I — k « I. 
Replacing the sum with an integral, using e~ x ~ (1 — x) 
and ( fc ) ~ l k , we obtain the approximate expression 

\C k \ oc N ini e k fl^e^dl (2) 

Jk 



When evaluating this integral, we have to consider 
three possible cases: 

1. The integral is independent of the cutoff because 
k < 7 — 1. In this case we obtain 

\C k \ ~N ini e k . (3) 

2. k > 7— 1 and e" 1 < fc ma x: m this case the exponen- 
tial function determines the cutoff to the integral, 
and we obtain 

\C k \ -AW 7 - 1 - (4) 

3. k > 7 — 1 and e _1 > k maK : in this case fc ma x deter- 
mines the cutoff to the integral, and we obtain 

|C fc |~A ini£ fc fc^ +1 . (5) 

The stochastic process ends when no nodes are left 
in container Co- On an average, the number of nodes 
in container Co is identical to the number of nonfrozen 
inputs minus the number of nonfrozen nodes, since the 
network is critical. If we neglect stochastic fluctuations 
during the process, the number of nodes in container Co 
becomes zero at the same time when the number of nodes 
in container Ck with k > 1 becomes zero, i.e. when 
e = 0. However, stochastic fluctuations will terminate 
the process earlier, at the moment where the fluctuations 
of the number of frozen nodes become of the same order 
as the expected number of frozen nodes. The variance of 
the number of frozen nodes is evaluated as follows: The 
probability that a given input has not yet become frozen 
at the moment where A nodes are left in the system, is 
e. When e is small, the number of nonfrozen inputs is 
Poisson distributed, with the variance being identical to 
the mean, which is proportional to A^e. For small e, the 
vast majority of nonfrozen inputs is found in container 
C\. Now a node in container C\ would be in container 
Co had its remaining input also become frozen during the 
process, and it follows that the variance of the number 
of frozen nodes is also of the order A; n ;e. The typical 
fluctuations in the number of frozen nodes are therefore 
of the order = n/A. Equating this number with 

the expected number of nodes in Co, which in turn is 
of the same order as the expected number of nodes in 
C2 , we obtain the following condition for the end of the 
stochastic process, where A is identical to the number of 
nonfrozen nodes, A n f: 

|C 2 | ~ = v/fe. (6) 

Depending on the value of 7 and on the dependence of 
fcmax on Ai n i, the number of nonfrozen nodes scales in a 
different way with Ai n j. 

For 7 > 3, the first of the three above cases applies to 
I C2 1 , and solving condition © for A n f we obtain 

Nnf ~ Af n f (7) 
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FIG. 2: Scaling collapse for the total number of nonfrozen nodes (x = 1) and for the number of nodes with two nonfrozen 
inputs (x — 2, left curves in each graph, shifted to the left by a factor 10 for better visibility), for three different values of 7 
and for the two different ways of choosing the input distributions. The function a(l,7,fc max ) is the appropriate exponent in 
Eqs. (7) to (9), and a(2,7,fc max ) = a(l, 7, fc max )/2. Each data set is generated by averaging over 10 5 realizations. 



at the end of the stochastic process. This is the same 
result as for a RBN with fixed k. Whenever the input 
distribution P(k) has a finite second moment, the num- 



ber of nonfrozen nodes scales as N- 



2/3 



and the number 



of nonfrozen nodes with two nonfrozen inputs scales as 



N, 



1/3 



The number of nonfrozen nodes with more than 



two nonfrozen inputs depends on whether k < 7 — 1, but 
it is in any case much smaller than the the number of 



nonfrozen nodes with two nonfrozen inputs, and we do 
not evaluate it here further. 

When 2 < 7 < 3 and when /c max oc iV in i, the second 
case applies, and we obtain using Eq. (|4|) 



(2 7 -4)/(2 7 -3) 



N nf ~ N[_ 

For 7 = 3, the exponent is 2/3, 
decreases to as 7 approaches 2. 



(8) 
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When 2 < 7 < 3 and when fc max oc (which is the 

case when the input distribution is fixed), the third case 
applies, and we obtain 



(9) 



These results Eqs. <j6j) - ((9j) should be also valid when 
biased Boolean functions are chosen. In this case, there 
is a nonvanishing probability that a node in container 
Cfe with k > 1 becomes frozen when an input becomes 
frozen. Therefore the expression {5} for \Ck\ obtains an 
additional factor 1 — p 2k — (1 — p) 2 , which is never close 
to and therefore does not change the scaling behavior 
of the integral. 

Our computer simulations confirm all these analytical 
considerations. As an example, we show in FigJ5] results 
obtained for the case of only constant and reversible func- 
tions, for both ways of choosing the input distributions. 
The excellent quality of the data collapses confirms our 
analytical calculations. 

Our results have a variety of implications. First, they 
show that many properties obtained for critical networks 
with a fixed number of inputs apply also to networks with 
a scale-free in-degree distribution once the frozen nodes 
have been removed. In particular, the number of non- 
frozen nodes with more than one nonfrozen input scales 
as the square root of the number of nonfrozen nodes. 
Only the dependence of the number of nonfrozen nodes 
on the total number of nodes is changed when 7 G (2,3 
We can therefore take over the results obtained in 
based on these properties of the nonfrozen nodes. It 
follows in particular that the number of relevant nodes 
in networks with a scale-free input distribution scales as 
the square root of the number of nonfrozen nodes, and 
that the number of relevant components is of the order of 
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log7Vi n i, with all but a limited number of relevant com- 
ponents being simple loops. It therefore follows again 
that the mean number and length of attractors diverges 
faster than any power law with the network size. This 
explains the finding in [l3T ] that the state-space structure 
of critical RBNs with fixed k and with a power-law input 
distribution is similar. Second, the number of nonfrozen 
nodes decreases with decreasing 7 £ (2,3), because the 
exponent becomes smaller. This explains why several 
authors have seen more frozen nodes and shorter attrac- 
tors in scale-free networks compared to standard RBNs. 
Third, the set of nonfrozen and relevant nodes is domi- 
nated by nodes with many inputs. This is due to the fact 
that each input has the same probability of surviving the 
stochastic process until the end. The average number 
of inputs of a node that has a surviving link is propor- 
tional to J k 2 N(k)dk, which is dominated by k max for 
7 € (2,3). When a relevant node is perturbed, the at- 
tractor is changed with a large probability, however when 
a frozen node is changed, the attractor changes with a 
probability that vanishes in the limit N — > 00. This ex- 
plains the findings in [12, HH, that attractors respond 
sensitively mainly to perturbations of highly connected 
nodes. Fourth, our results disagree with the finite-size 
arguments in [l8| , which predict that the number of non- 
frozen nodes scales as -N;^ • This is in our view due 
to the fact that an infinite (sustained) perturbation has 
properties that are fundamentally different from those of 
finite perturbations, in which case arguments based on 
finite-size scaling do not work. 
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